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Abstract 

We study a dispersive counterpart of the classical gas dynamics problem of the 
interaction of a shock wave with a counter-propagating simple rarefaction wave often 
referred to as the shock wave refraction. The refraction of a one-dimensional dispersive 
shock wave (DSW) due to its head-on collision with the centred rarefaction wave (RW) 
is considered in the framework of defocusing nonlinear Schrodinger (NLS) equation. 
For the integrable cubic nonlinearity case we present a full asymptotic description 
of the DSW refraction by constructing appropriate exact solutions of the Whitham 
modulation equations in Riemann invariants. For the NLS equation with saturable 
nonlinearity, whose modulation system does not possess Riemann invariants, we take 
advantage of the recently developed method for the DSW description in non-integrable 
dispersive systems to obtain main physical parameters of the DSW refraction. The key 
features of the DSW-RW interaction predicted by our modulation theory analysis are 
confirmed by direct numerical solutions of the full dispersive problem. 



1 Introduction 



Recent developments of experimental techniques of cold-atom and laser physics and ob- 
servations of a number of superfluid and optical counterparts of classical hydrodynamic 
phenomena such as solitons, shock waves, rarefaction waves, vortex streets etc. (see e.g. 
[SI [311 |55l [22l [191 Eni [3]) stimulated the growing interest in the mathematical methods and 
results of dispersive hydrodynamics — the theory of multiscale nonlinear flows in media with 
dispersive (rather than dissipative) mechanisms of regularization of breaking singularities. 
Central to dispersive hydrodynamics is the theory of dispersive shock waves (DSWs), which 
represent expanding nonlinear wavetrains connecting two different hydrodynamic states and 
replacing, in conservative continuous media, the classical viscous shocks (see [30] and refer- 
ences therein). Owing to their rich dynamics and fundamental physical nature, the DSWs 
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have recently become an object of very active theoretical and experimental investigations, 
most notably in Bose-Einstein condensates (BECs) (see, e.g., [HI Ell E]), where these waves 
represent a striking manifestation of quantum statistics on a macroscopic scale. Another 
area of active modern DSW research is nonlinear optics (see [391 EH l22l |5]). 

While the dynamics of isolated DSWs have been studied in numerous works since the 
pioneering paper [27j by Gurevich and Pitaevskii, their interaction behaviour has begun to 
be investigated relatively recently. One can distinguish two prototypical problems arising 
in this connection: interaction two DSWs and interaction of a DSW with a simple (rar- 
efaction) wave. Both problems admit full analytical description in the framework of the 
Whitham modulation theory associated with integrable wave equations, however, the former 
problem involves complicated analysis of nonlinear multiphase wavetrains (see e.g. [23] for 
the KdV equation) so in applied problems one is usually better off finding the corresponding 
modulation solutions numerically (see [1], [6], [29j). The latter problem of the DSW-RW 
interaction, on the contrary, involves only single-phase dynamics, so one could hope for a 
relatively simple effective analytic description of the interaction. 

In [1] a complete classification of unidirectional interactions of DSWs and RWs in weakly 
dispersive flows was made using the analytical inverse scattering transform (1ST) solutions 
for the KdV equation and numerical solutions of the KdV- Whitham equations. This classi- 
fication has revealed certain similarities as well as fundamental differences between classical 
and dispersive-hydrodynamic overtaking shock wave-rarefaction wave interactions (see also 
[T5] for the analytical modulation solutions describing some of the KdV overtaking DSW- 
RW interactions). In many physical settings, however, one has to deal with bi-directional 
(head-on) wave collisions which cannot be captured by the KdV type models and should be 
studied in the framework of appropriate two-wave equations. Such a bi-directional DSW-RW 
interaction represents a dispersive counterpart of the classical gas dynamics problem which 
is often referred to as the "shock wave refraction" . 

When a one-dimensional viscous shock wave (SW) undergoes a head-on collision with a 
rarefaction wave, the parameters of two waves alter so that the long-time output of such an 
interaction consists of a new pair of SW and RW propagating in opposite directions. Since 
the SW speed changes from one constant value to another as a result of its propagation 
through the finite RW region with varying density and velocity, the interaction diagram in 
the (xt)-plane could be naturally interpreted as the SW refraction on the RW. As a matter 
of fact, the SW refraction can be observed in two-dimensional stationary flows where the 
effect acquires its direct geometrical significance as deflection of the SW from its original 
propagation direction accompanied by the change of its strength and speed. 

Refraction of SW's has been the subject of many gas and fluid dynamics investigations 
(see, e.g. the original wartime report [10] by Courant and Friedrichs and some of the well- 
known research papers [lH IH |28l HH] as well as classical monographs [11], [19], [15]). It must 
be said, however, that, while the qualitative features of the SW refraction process have been 
understood very well, its analytical description is seriously hindered due to the presence of 
the varying entropy region between the refracted SW and RW. As a result, the system of 
equations describing the head-on SW-RW interaction turns out to be so complicated that 
numerical solution becomes in most cases the only available resort. 

In dispersive compressible dissipationless flows the entropy does not change and, in con- 
trast to viscous gas dynamics, the bidirectional DSW-RW interaction can be described an- 
alytically in terms of solutions of the Whitham modulation equations [58j associated with 
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the original dispersive-hydrodynamic system and governing slow variations of the wave pa- 
rameters (amplitude, wavenumber, mean etc.) on the scale much larger than the medium 
typical coherence length. 

In this paper, we perform an analytical study of the head-on DSW-RW interaction in 
the framework of the nonlinear Schrodinger equation with defocussing, which is a standard 
mathematical model in nonlinear optics and condensed matter physics (see e.g. [38], |47j ) . 
Thus, apart from the obvious theoretical significance as a dispersive counterpart of a classical 
gas dynamics problem, the theory of the DSW refraction in nonlinear Schrodinger flows is 
fundamental to the understanding of 'dispersive-hydrodynamic' flow interactions in super- 
fluids and nonlinear optical media. For the case of cubic nonlinearity the NLS equation is a 
completely integrable system and a full asymptotic description of the DSW-RW interaction 
becomes possible owing to the availability of exact solutions of the NLS-Whitham equations 
describing slow variations of the rapidly oscillating wave field in the interaction zone. The 
key element of the analytical construction is the mapping of the two-component reduction 
of the NLS-Whitham system to the classical linear Euler-Poisson-Darboux (EPD) equation. 
This mapping was introduced for the KdV-Whitham system in [Jl] [25], [26] and [51]; and for 
the NLS equation in [26]. Remarkably, the same EPD equation describes, on the hodograph 
plane, the interaction of two nonlinear simple waves in ideal shallow-water dynamics - see, 

e.g. m- 

Along with the study of the DSW refraction in Kerr media described by the integrable 
NLS equation, we also undertake a similar investigation of the DSW-RW interaction in 
the framework of the NLS equation with saturable nonlinearity (sNLS), which represents a 
standard model for the optical beam propagation in photorefractive crystals (see, e.g. [2Tj , 
[8], [3H]). The photorefractive systems have been recently used for the modelling dispersive- 
hydrodynamic flows in BECs by means of an all-optical setting [55] so the quantification of 
the contribution of the saturation effects to the 'superfluid' dynamics of light is important 
for the comparison with BEC experiments. 

The sNLS equation is not integrable by the inverse scattering transform, so the associated 
Whitham system does not possess Riemann invariants. As a result, this system cannot be 
reduced to the EPD equation and the analytic method employed for the description of the 
DSW refraction in the cubic nonlinearity case is not applicable to the sNLS equation. To 
tackle the sNLS refraction problem analytically, we take advantage of the approach to the 
dispersive Riemann problem treatment in non-integrable conservative systems developed in 
[12J . This has enabled us to derive analytically the key parameters of the refracted DSW, 
as well as the DSW refraction, acceleration and amplification coefficients as functions of 
the initial data and the saturation parameter 7. We note that the theory of propagation 
of simple photorefractive DSWs was developed in paper [13], which contains some detailed 
explanations of the application of the method of [12] to the sNLS equation. In the present 
paper, we extend the results of [13] to describe the photorefractive DSW-RW interaction. 
In particular, we show that for a broad range of parameters the photorefractive DSW-RW 
interaction is asymptotically "clean", i.e. is not accompanied by the generation of new DSWs 
or RWs. 

The direct numerical simulations for NLS and sNLS equations (using standard split-step 
Fourier method - see, e.g., pTj) confirm all the key features of the bidirectional DSW- 
RW interaction, predicted by our modulation analysis. We stress, however, that, while we 
perform some basic comparisons for the typical behaviours of the key physical parameters, a 
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systematic validation of the obtained solutions in the framework of full dispersive problem is 
beyond the scope of the present paper and would require, in particular, a detailed comparison 
with the numerical solutions of the small-dispersion limit of the NLS equation (see, e.g., [SS] 
for the corresponding analysis for the KdV equation). At the same time we would like to 
mention that there have been a number of comparisons between the solutions of the NLS- 
Whitham systems and direct numerics for the NLS equation in recent literature (see e.g. 
[36] . [3T] . |13] . [T7]). all of them showing a very good agreement, so we have some confidence 
in the relevance of our modulation solutions to the properties of the full dispersive DSW-RW 
interaction dynamics. 



2 DSW refraction in Kerr media: formulation of the 
problem 

We first formulate the problem for the defocusing NLS equation with cubic (Kerr) nonlin- 
earity 

ieiJt + -^^^-\^\^ij = 0, (1) 

where ip is a complex valued function and e is a dimensionless dispersion parameter (coher- 
ence length). Using the Madelung transformation ip i— )■ {n,u) 



ipix, t) = \/ n{x, t) exp ( - / u{x',t)dx'\ , (2) 



e 

where n{x,t) > and u{x,t) are real-valued functions, we represent the NLS equation ([1]) 
in the "dispersive-hydrodynamic" form 

rit + {nu)^ = 0, 
2 I nl nxx\ „ (3) 



ut + uu^ + n^ + e — - -;— =0 

with the 'fluid' density n and velocity u. 

The dispersionless (classical) limit of system ([3]) is obtained by setting e = and is 
nothing but the system of ideal shallow-water equations 

fit + {nu)x = 0, ut + uux = , (4) 

which can be represented in the diagonal form 

^ + t/,(A,.A_)— = 0, (5) 

with the Riemann invariants 



\± = ]^u±^/n (6) 



and the characteristic velocities 
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To study the bidirectional (head-on) interaction of a DSW and a rarefaction wave (RW) 
we consider the following configuration. Let at some moment of time say t = tc, & simple 
right-propagating DSW confined to the expanding region x^it) < x < xf(t) and a simple 
left-propagating RW located at X2(t) < x < be separated by an undisturbed flow 

region x^it) < x < X2{t) with n = 1 and u = (see Fig. 1). Without much loss of 
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Figure 1: Sketch of the density profile in the NLS flow prior to head-on DSW-RW interaction 



generality one can assume that the DSW and RW are both centred in the (x, t)-plane at 
(0, 0) and (0, /) respectively, so that xf = sft, xf = l + sft, where s]^ > > 0, < < 
are the speeds of the respective DSW and RW edges. We also assume that / ^ 1. 

The following transition conditions must be satisfied across the DSW and RW respectively 

(see 



A+(X2 , tc) 



\^{x^,tc) 
A_(,T2 



^2 5 ^cj 



simple right-propagating DSW transition (8) 
simple left-propagating RW transition (9) 



The transition conditions (jH]), (E]) imply that the described above flow configuration can be 
realised as a result of the evolution of the initial flow profile n{x, 0), u{x, 0) specified in terms 
of the shallow- water Riemann invariants X± (Q having the jumps of different polarity shifted 
with respect to one another by a distance / (see Fig 2a): 
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Figure 2: Initial conditions for the NLS equation ([3]) leading to the head-on DSW-RW 
interaction. Left: hydrodynamic Riemann invariants \± ( !T0|) : Right: corresponding density 
n and velocity u distributions ( ITTj) . 
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\ t n\ \ for X < 0, \ / n^ / X<1, , . 

^+^^'0) = \l for x>0; A_(x,0) = |^_ for x > ^^^^ 

where ^4+ > 1 and -I < < 1. 

The initial conditions for n and m corresponding to ffTOj) are then readily found using ([6]) 
in the form of piecewise constant distributions (see Fig. 2b) 

( \{l + A+f>l for a;<0, f v4+ - 1 > for x < 

n(x, 0) = < 1 for < s < /, u{x, 0) = < for < x 

\\{l-A-f<l ioi x>l- [l + A->0 for x>/ 

(11) 

Our concern will be to obtain analytical description of the head-on DSW-RW interaction 
in terms of the initial profile parameters A"*", A~ and /. 

The evolution ([3]), ( ITTj) can be qualitatively understood using the results of papers [H] 
and [23] where the Riemann problem (which is a particular case of the problem (|3]), ( !T0|) 
with / = 0) was considered and a full classification for the different cases of the decay was 
constructed using similarity solutions of the modulation NLS-Whitham equations in the 
framework of the "matched regularisation" procedure of the Gurevich-Pitaevskii type (see 
also [S] for the further detailed analysis using the alternative "global regularisation" 
formulation). The crucial difference between the dispersive Riemann problem of |T3], |24] . 
and the present problem ([3]), f lTTj) is that the two discontinuities for A+ and A_ are now 
spaced by a large distance / so the modulation problem is no longer self-similar and a more 
general consideration is required. 

The asymptotic solution of the NLS dispersive Riemann problem obtained in [21], [ll] 
(see also [39], [6]) and our direct numerical simulations of the more general general initial- 
value problem ([3]), ( [TT]) suggest that the evolution (|3]), ( ITTj) will initially lead to the formation 
of a right-propagating simple DSW and a left-propagating simple RW as in Figs. 1,2. Both 
waves expand with time and start to overlap and interact at some t = t^. The interaction 
continues until some t = t* > Iq when the two waves fully separate so that at t > t* there 
is a combination of new, "refracted", simple DSW and RW separated by a new constant 
state rio 7^ 1, tto 7^ 0. All the described stages of the DSW refraction are clearly seen on 
the direct numerical simulation plots in Fig. 3, 4. The most obvious effect of the head-on 
DSW-RW interaction seen in Fig. 3 is the change of the key parameters (intensities, speeds) 
of the interacting waves. Another, more subtle, effect is the change of the phase distribution 
acquired by the DSW during the interaction. The refraction phase shift d of the DSW trailing 
dark soliton is shown in Fig. 4. Also we note that, if the incident DSW (RW) was centred 
at if: = 0, the refracted DSW (RW) will generally no longer be a centred wave. In Section 5 
we shall construct an exact analytic solution of the NLS-Whitham equations asymptotically 
describing all stages of the head-on DSW-RW interaction seen in Figs. 3, 4 and compare the 
definitive DSW refraction parameters (the DSW amplification and acceleration coefficients as 
well as the refraction phase shift) with the corresponding parameters obtained numerically. 

In conclusion of this section we note that the outlined head-on DSW-RW interaction can 
be naturally realised in the framework of the dispersive piston problem (see [32], [17], [37] ) 
involving two pistons, the right piston being pulled out from the gas with constant velocity 
producing thus a left-propagating rarefaction wave while the left piston being pushed into 
the gas producing the right propagating DSW. Another pertinent problem is the interaction 
of stationary two-dimensional DSW and RW forming in hypersonic dispersive flows past 
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Figure 3: Bidirectional interaction of a DSW and RW: density (upper) and velocity (lower) 
profile; Initial data parameters: = 1.5, A~ = —0.4, / = 50. The value of the dispersion 
parameter e used in the simulations is 0.4 
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Figure 4: (Colour online) Density plot corresponding to the DSW-RW interaction shown in 
Fig. 3. The regions are as follows: I - incident DSW; II - incident RW; III - DSW-RW 
interaction region; IV - Refracted RW; V - Refracted DSW. 

extended obstacles. This latter configuration is relevant to the BEC experiments j9j and can 
also be reformulated in terms of the already mentioned dispersive piston problem (see jl6] . 

m)- 

3 Refraction of shock waves in classical gas dynamics 

Before we proceed with the analysis of the bidirectional dispersive refraction problem (jS]), ( ITTl) 
we outline some classical results on the head-on interaction of viscous shocks and rarefaction 
waves (see e.g. [10], gl], |19], |15]). 

Consider a one- dimensional motion of a polytropic isentropic gas, i.e. a gas with the 
equation of state p = crC , where p and n are the gas pressure and density respectively, 7 is 
the adiabatic exponent and c is a constant (the dispersionless shallow-water dynamics (jlj) is 
equivalent to the dynamics of the polytropic gas with 7 = 2). We consider the following flow 
configuration (see Fig. 5). Let the gas motion at some moment of time, say t = tc > consist 
of three regions of constant flow separated by two waves: a right-propagating shock wave 
(SW) located at some x = Xc and a left-propagating RW centred aX x = I and occupying 
a finite region of space (as already was mentioned, such a configuration can be created by 
piston motion inside a tube - see e.g. Let the density and velocity of the flow be 

{rii.ui) as a; — 7- —00 and {n2,U2) as x ^ -|-oo. Then the gas motion at t > tc can be 
qualitatively described as follows: 

• The SW and RW propagate independently until the moment t = to, when the shock 
enters the rarefaction wave region at some x = Xq say. Before that moment, i.e. for < 
t < to, the entropy undergoes a rapid constant change across the SW so the SW speed 
and strength (the pressure excess across it) are determined by the standard Rankine- 
Hugoniot conditions. The RW is described by the centred left-propagating simple- 
wave solution of the inviscid hydrodynamic equations of motion. The parameters of 
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Figure 5: Head-on interaction of SW and RW in classical gas dynamics 



the constant flow between the SW and RW are found at the intersection of the n-u 
diagrams for the SW and RW (see e.g. fl]). 

• During certain time interval to < t < t* the SW and the RW interact. The interaction 
is accompanied by the variations of the shock strength and results in the formation of 
the varying entropy region (the so-called 'entropy wave') behind the SW. Therefore, 
the flow behind the refracted SW is not isentropic. 

• At t = t* the SW exits the RW region and the two waves again propagate separately 
in opposite directions, each having an altered (as compared with the values before the 
interaction) set of parameters. An important general result is that the speeds of the 
refracted SW and RW and the density /velocity jumps across them are exactly as they 
would have been in the corresponding origin-centred Riemann problem (i.e. in the 
decay of an initial discontinuity problem with the gas parameters (ni, ui) at x < and 
(^2,^2) at X > ), however, the spatial locations of the refracted waves differ from 
those in the corresponding Riemann problem. The refracted SW always has greater 
speed and strength than the original one. 

As already was mentioned, the presence of the 'entropy wave' behind the refracted SW rad- 
ically complicates quantitative analysis of the motion and, as a result, the SW-RW head-on 
collision problem can generally be treated only numerically (we stress that the notion of the 
entropy wave applies to viscous SWs in gas dynamics and is not applicable to, say, shal- 
low water dynamics). In contrast to classical gas dynamics, dispersive- hydrodynamic flows 
governed by completely integrable equations often admit full analytical description. In par- 
ticular, such a description is available for the DSW refraction process. This description can 
also be generalised (to some extent) to certain types of non-integrable dispersive equations. 

4 Single-phase modulation theory for the defocusing 
cubic NLS equation: account of results 

It is known very well that analytical theory of one- dimensional dispersive compressible flows 
containing DSWs can be constructed in the framework of the Whitham modulation equations 
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[58] . In this section we make a brief account of the relevant results of the modulation theory 
for the defocusing cubic NLS equation, which will be necessary for the analysis of the DSW 
- RW interaction in the subsequent sections. The single-phase NLS modulation system 
was derived in [20], [16] (see also [39]) using the finite-gap integration methods. A more 
elementary derivation of this system using a reduced version of the single-gap integration 
can be found in [35]. Importantly, the theory presented in this section makes substantial 
use of the integrability of the NLS-Whitham modulation system, which is inherited from 
the complete integrability of the original cubic NLS equation. A different method, proposed 
in [12] and applicable to the description of DSWs in nonintegrable systems, will be used 
in Section 6 for the description of the DSW refraction in the media described by the NLS 
equation with saturable nonlinearity (!97|) . which does not enjoy the complete integrability 
property. 

It should be noted that, since the results of the modulation theory do not depend on 
the value of the dispersion parameter e in the NLS equation ([T]), we shall assume e = 1 in 
the subsequent analytical representations of the periodic solutions, while in the numerical 
simulations we shall normally be using smaller values of e to reduce the temporal scale of 
the DSW structure establishment. 



4.1 Periodic solution and modulation equations 

The periodic travelling wave solution of the defocusing NLS equation ([3]) can be expressed 
in terms of the Jacobi elliptic sn function and is parametrised by four integrals of motion 
Ai < A2 < A3 < A4 (as was already mentioned, we assume e = 1 in the NLS equation), 

n = i(A4 - A3 - A2 + Ai)2 + (A4 - A3)(A2 - Ai) sn2 (v/(A4 - A2)(A3 - Ai) 0, m) , (12) 

u = U--, (13) 
n 

where C = i(-Ai - A2 + A3 + A4)(-Ai + A2 - A3 + A4)(Ai - A2 - A3 + A4), 

e = x-ut-eo, u = -J2^^^ (14) 

i=l 

U being the phase velocity of the nonlinear wave and 6q the initial phase. 
The modulus < m < 1 of the elliptic solution f|T2l) is defined as 

(A2-A0(A4-A3) 

(A4-A2)(A3-A0' ^ ^ 

and the wave amplitude is 

a = (A4-A3)(A2-Ai). (16) 
The wavelength of the periodic wave (fT2|) is given by 



A4 A2 

dx r dx 



. v/(A-Ai)(A-A2)(A-A3)(A4-A) J ^(A - Ai)(A2 - A)(A3 - A)(A4 - A) 

A3 Ai \^') 

2K(m) 



V(A4-A2)(A3-Ai)' 
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K(m) being the complete elhptic integral of the first kind. As a matter of fact, £ > 0. 

In the limit as m — )■ 1 (i.e. as A3 — )■ A2) the travelling wave solution (IT^ turns into a 
dark soliton ^ 

n = ns- ,2 / ^ , (18) 



cosh {y/a^{x — Ugt — Oq)) 

where the background density n^, the soliton amplitude and velocity Ug are expressed in 
terms of Ai, A2, A4 as 

ris = ^(A4 - as = (A4 - A2)(A2 - Ai), = ^(Ai + 2A2 + A4) . (19) 

Allowing the parameters Ai, A2, A3, A4 of the travelling wave solution ( !T2|) to be slowly vary- 
ing functions of x and t, one arrives, via the averaging or an equivalent multiple-scale per- 
turbation procedure, at a modulated nonlinear periodic wave in which the evolution of 
A = {Ai, A2, A3, A4} is governed by the Whitham modulation equations [201 SS] (see [SS, ES] 
for a detailed description of the Whitham method) 

^ + ^.(A)^ = 0, ^ = 1,2,3,4, (20) 

Aj's being the Riemann invariants. The characteristic velocities can be computed using the 
formula [2Sl ES] 

l^i(A) = (^1 - ^a,^ f/, z = 1,2,3,4, where d, = d/d\. (21) 
Substitution of Eq. ( |T7I) into Eq. ( 12T|) yields the explicit expressions 







V2 = 




V3 = 











(A4- 


Ai)(A2 


-Ai) 


(A4 


-Ai) 


-(A4- 


- A2)Ai(m) 




(A3- 


A2)(A2 


-Ai) 


(A3 


-A2) 


-(A3- 


- Ai)/i(m) 




(A4- 


A3)(A3 


-A2) 


(A3 


-A2) 


-(A4- 


- A2)Ai(m) 




(A4- 


A3)(A4 


-Ai) 


(A4 


-Ai) 


-(A3- 


- Ai)/i(m) 



(22) 



where yu(m) = E(m)/K(m), E(m) being the complete elliptic integral of the second kind. 
The characteristic velocities ( 122|) are real for all values of the Riemann invariants, therefore 
system (120|) is hyperbolic. Moreover, it is not difficult to show using representation (12T]) that 

diVi > for all i , (23) 

so the NLS- Whitham system (l20l) . fl22|) is genuinely nonlinear [53]. Indeed, differentiating 
dni) we get: 

a.v. = ^au. (24) 

Using the integral representations (fTTj) for £ one can readily see that df^Q > for all i (it is 
convenient to use the first representation for the differentiations with respect Ai and A2 and 
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the second one for the differentiations with respect A3 and A4), which immediately imphes 



Also, using f l21l) and the intergral representations (IT7|) one can readily show by a direct 
calculation that 

i > j implies Vi > Vj . (25) 

Thus, the ordering Ai < A2 < A3 < A4 of the Riemann invariants implies a similar ordering 
^1 < ^2 < V3 < V4 for the characteristic velocities. We note that the properties (l23ll and (!25|) 
were established in [31], [39] using the finite-gap integration framework for the derivation of 
the Whitham equations. 

For the DSW analysis in the subsequent sections we shall need the reductions of formulae 
( I22I) for the limiting cases m = (harmonic limit) and m = 1 (soliton limit). 

The harmonic limit m = can be achieved in one of the two possible ways: either via 
A2 = Ai or via A3 = A4. Then: 



K X X T/ T/ X , A3 + A4 , 2(A3-Ai)(A4-Ai) 

when A2 = Ai : Vo = Vi = Ai H ; ; ; 

2 2A1-A3-A4 



^3 = ^A3 + ^A4 = V-(A3, A4) , = ^A4 + ^A3 = V+iXs, A4) 

>, X X T/ T/ X , A1 + A2 , 2(A4-A2)(A4-Ai) 
when A3 = A4 : V3 = V4 = A4 H h 



(26) 



2A4 — A2 — Ai 



(27) 



V, = ^Ai + ^A2 = V^Xi, A2) , V2 = + ^Ai = r+(Ai, A2) . 
In the soliton limit we have m = 1. This can happen only if A2 = A3, so we obtain: 

when A2 = A3 : V2 = V3 = ^(Ai + 2A2 + A4) = U, , 

Vi = ^Ai + ^A4 = K.(Ai, A4) , V, = ^A4 + ^Ai = y+(Ai, A4) . 



(2^ 



Thus, in both harmonic and soliton limits the fourth-order modulation system fl20|) . f l22|) 
reduces to the system of three equations, two of which are decoupled. Moreover, one can see 
that in all considered limiting cases the decoupled equations agree with the dispersionless 
limit of the NLS equation ([1]). This property makes possible the matching of the modulation 
solution with the solution to the dispersionless limit equations at the points where m = or 
m = 1. 



4.2 Free-boundary matching conditions for the modulation equa- 
tions 

In the description of a DSW, the Whitham equations ( 120|) must be endowed with certain 
initial or boundary conditions for the Riemann invariants Aj. We shall be using the Gurevich- 
Pitaevskii type boundary-value (matching) problem first formulated in [27] for the KdV dis- 
persive shock waves and extended to the NLS case in [21]. A different type of the problem 
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formulation (the so-called regularised initial- value problem for the NLS-Whitham equations) 
proposed in [39] and recently used in [B], [53] for the numerical analysis of the DSW interac- 
tion is less convenient for our purposes as the analytical description of the interaction zone 
requires the hodograph solutions of the Whitham equations, and the poor compatibility of 
the initial- value problems with the hodograph transform is well known (see e.g. [58] )■ The 
Gurevich-Pitaevskii matching conditions, on the contrary, are ideally compatible with the 
hodograph transform as they turn into the classical Goursat type characteristic boundary 
conditions on the hodograph plane (see [25], [26], [18], [IZ]). It is clear that both formulations 
(regularised initial-value problem for the Whitham equations and the Gurevich-Pitaevskii 
type matching problem) must be equivalent, although we are not aware of the rigorous proof 
of this equivalence. 

To be specific, we shall formulate boundary (matching) conditions for the right-propagating 
DSW. Without loss of generality we assume that the formation of the DSW starts at the 
origin of the (x, t)-plane. In the Gurevich-Pitaevskii setting, the upper (x, t)-half plane is 
split into three regions (see Fig. 6): (— oo,x~(t)), [x^ {t) , {t)] and (x+(t),+oo). 




Figure 6: Splitting of the xt-plane in the Gurevich-Pitaevskii problem for the defocusing 
NLS equation. 

In the "outer" regions {—oo,x~(t)) and {x^{t), +oo) the flow is governed by the disper- 
sionless limit of the NLS equation, i.e. by the shallow- water system ([5]), for the Riemann 
invariants X±. In the DSW region [x~ (t) , x^ (t)] the averaged oscillatory flow is described 
by four Whitham equations (!20|) for the Riemann invariants Xj with the following matching 
conditions at the trailing x~(t) and leading x^{t) edges of the DSW (see [211 [18] for details): 

At X = x"(t) : As = A2 , A4 = A+, Ai = A_ , , . 

At x = x+(t) : A3 = A4 , A2 = A+, Ai = A_ . ^ ' 

Here X±{x,t) are the Riemann invariants of the dispersionless limit of the NLS equation in 
the hydrodynamic form (E]), ((Tj). The free boundaries x^{t) of the DSW are defined by the 
kinematic conditions 

dx~ dx^ 

= V2(Ai, A2, A2, A4) = V3(Ai, A2, A2, A4) , = V3(Ai, A2, A4, A4) = V4(Ai, A2, A4, A4) 

(30) 

and so are multiple characteristics of the Whitham system. The multiple characteristic veloc- 
ities V2 = V3 and V3 = V4 in ( 130|) are explicitly given by equations (128|) and (!27|) respectively. 
One should stress that determination of x^{t) is an inherent part of the construction of the 
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full modulation solution. We also emphasize that matching conditions (129|1 are consistent 
with the structure of the Whitham system f l20p . (12^ in the limiting cases m = and m = 1 
(see (125]) ) and with the spatial oscillatory profile of the DSW in the defocusing NLS 

dispersive hydrodynamics (as is known very well, such a DSW has a dark soliton (m = 1) 
at the trailing edge and degenerates into the vanishing amplitude harmonic wave (m = 0) 
at the leading edge— see [211 [12 [36l |3T]). 



4.3 Hodograph transform and the mapping to the Euler-Poisson- 
Darboux equation 

The hydrodynamic type modulation system (120|) . (122|) can be reduced to a system of linear 
partial differential equations using the (generalised) hodograph transform [51]. We first fix 
two of the Riemann invariants, say 

K = = constant , Xj = Xjo = constant , i ^ j , (31) 

to reduce fl20|) to the system of two equations for the remaining two invariants Xk{x,t) and 

Xi{x,t), kj^ly^iy^j 

where Vk,i{Xk, Xi) = Vkj{Xio, Xjq, Xk, Xi). Applying the hodograph transform to system fl32|) 
we arrive at a linear system for x{Xk, Xi), t{Xk, Xi), 

Vi{X,,Xi)^ = 0, |^-r,(A,,A0# = O. (33) 



oXk oXk oXi dXi 

Note that the hodograph transform requires that d^Xk^i 7^ 0. Now we make in fl55]) the 
change of variables 

x-Vut = Wk, x-Vit = Wi, (34) 
which reduces it to a symmetric system for Wk{Xk, A/), Wi{Xk, Xi): 

k^l; dk = d/dXk . (35) 



Wk -Wi Vk- Vi 

The symmetry between Vi and Wi in ( l35l) and the 'potential' structure ( 12T]) of the vector 
function (V^, VJ) implies the possibility of introducing a single scalar function g{Xk, Xi) instead 
of the vector iWk,Wi): 

= {^'-^^^3, i = kj, (36) 

or, which is the same (use (121])). 

W, = g + 2{V,-U)d,g, i = kj . (37) 



Then substituting f l2T|) . fl36|l into f l35|) we arrive, taking into account ( !T7|) . at the Euler- 
Poisson-Darboux (EPD) equation for g{Xk, A;), 

2(A;-Afc)52,^ = %-9fc^. (38) 
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The EPD equation was first derived in the present NLS context in [2^ and was later used 

in [TS] and [22] for the construction of the general solution of the semi-classical defocusing 
NLS equation with smooth monotonically decreasing initial data. 

Note that system fl32|) essentially describes the interaction of two simple waves of modu- 
lation of the NLS equation so the fact that this system reduces, in the hodograph plane, to 
the EPD equation fl38l) describing interaction of two simple waves in classical dispersionless 
shallow-water theory (or in gas dynamics of polytropic isentropic gas with 7 = 2) (see [58] 
for instance) is quite remarkable. 

The general solution of the EPD equation (138|1 can be represented in the form (see, for 
instance, [53] ) 



where 0i,2(A) are arbitrary (generally, complex- valued) functions and ai_2 are arbitrary con- 
stants (which could be absorbed into (f>i,2)- 

As a matter of fact, the same construction can be realized for any pair of Riemann 
invariants while the two remaining invariants are fixed. Moreover, equations f p4|) - fl35|) and 
further fl^T]) - fl55]) turn out to be valid even when all four Riemann invariants vary [2B| [T5]. 
This becomes possible for two reasons. Firstly, the NLS modulation system fl2U]) . (12T]) is 
integrable via the generalized hodograph transform [51] which reduces it to overdetermined 
consistent system ( 135|1 . where k,l = 1,2,3,4, k ^ I. Secondly, the "potential" structure 
of the characteristic speeds ( 121]) makes it possible to use the same substitution (!36|) for all 
= 1,2,3,4 which results in the consistent system of six EPD equations f l38p involving all 
pairs Afc, A;, k I. 

Thus, the problem of integration of the nonlinear Whitham system fl20|) with rather com- 
plicated coefficients fl22|) is essentially reduced to solving the classical linear EPD equation 
fl38|) . Essentially, one needs to express the functions 0i,2(A) in the general solution f p9|) 
in terms of the initial or boundary conditions for the NLS equation ([T]). As was shown 
in [2S], [IB] (see also [T7]) the free-boundary nonlinear matching conditions fl29]) are most 
conveniently translated into a classical linear Goursat characteristic boundary problem for 
the EPD equation ( 138|) . This enables one to find the unknown functions 0i,2(A) in terms of 
Abel transforms of the initial data. 

In conclusion of this section we note that hodograph solutions do not include the special 
family of the simple-wave solutions as the latter correspond to the vanishing Jacobian of the 
hodograph transform (A^, A;) 1— {x,t) (see, for instance, [5S]). 

4.4 Modulation phase shift 

In the modulated wave, the initial phase ^0 of the periodic solution f[T^ - f lT5]) is no longer 
an independent constant parameter but rather a slow function of x, t so it is better described 
as the modulation phase shift. As was shown in [ID], the function 6o{x, t) can be found from 
the requirement that the local wavenumber k = 27r/£ and the local frequency u = kU in 
the modulated wave (fT2|) must satisfy the generalised phase relationships 




(39) 



k = Q 



uj = -e 



(40) 
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where 

Q = k9 = kx-ujt~ k9o (41) 
is the angular phase. Relationships fl4Up imply the 'conservation of waves' law 

kt + u, = 0, (42) 

which is consistent with the modulation system ( l20l) and thus yields the representation 
Vi = diu/dik for the characteristic speeds, equivalent to fl^ . 

For the general modulation relationship (HOj) to be consistent with the linear x, t-dependence 
of the phase pij) entering the local single-phase NLS solutions (!T2|) . (fT5|) one must assume 
9o{x,t) = '(9o(Ai, A2, A3, A4), which implies that the phase shift is completely determined 
by the evolution of the Riemann invariants Xj{x,t) in the modulation solution. To find 
''^o(Ai, A2, A3, A4) we differentiate (HTl) with respect to x to obtain 

4 

= k + ^ {xdik — tdiU — ^odik — kdi^o} dxK ■ (43) 

i=l 

Comparing ( l43l) with ( HOl) we obtain for any pair i,j, i ^ j 

xdik — tdiuo — "dodik — kdi^Q = , xdjk — tdju — d^djk — kdjdo = , (44) 

provided d^Xij 7^ 0. On using diu/dik = Vi and k = 2tx / 2, system (jH]) is readily transformed 
to the form 



x-Vnt= [I- ^-^dn]^Q, n = i,j, i^j. (45) 



Comparison of expression fj45l) with the modulation hodograph solution fl34l) . fl36l) enables 
one to identify the modulation phase shift ^o{X) = 9o{x,t) with the solution g{X) to the 
relevant boundary value problem for the EPD equation (|38|) . i.e. 



9o{x,t)=g{X{x,t)). (46) 

One can also see from P5|) that one should set ^0 = for a simple centred DSW described 
by the modulation solution in which all but one Riemann invariants are constants and the 
varying invariant, say A^, is implicitly specified by the equation x — Vmt = 0. The condition 
^0 = then implies that in the dispersive Riemann (decay of a step) problem the DSW 
trailing dark soliton f|T8l) is centred exactly at the trailing edge x~{t) defined by fl30|) . fl28l) . 

5 Interaction of DSW and RW: modulation solution 

5.1 Before interaction, < t < to 

At t = 0, a simple origin-centred right-propagating DSW is generated due to the jump of the 
Riemann invariant A+ while the jump of A_ produces a similarity "shallow-water" rarefaction 
wave centred at x = / and propagating to the left. 

The similarity modulation solution describing the DSW has the form [21], [H] 

Ai = -1, A2 = l, A4 = A+, 

^-I/r 1 1 A 4+l-^i±^ (-4^-A3)(A3-l) (47) 
j-Vs{-lA,X3,A )-^ A3 - 1 - (A+ - lUm) ' 
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Figure 7: Schematic behaviour for the Riemann invariants before the interaction of the 
DSW and RW, < t < to- 



where 

m = - (48) 

The boundaries of the DSW are then found from (147|) by setting A3 = (i.e. m = 0) for 
the leading edge x\ and A3 = 1 (i.e. m = 1) for the traihng edge x^: 

xr = ^^t, xt = {2A^-^)t. (49) 

The dark sohton at the traihng edge of the DSW has the amphtude and rides on the 
background rig defined by (see ( |T9l) ) 



= (1 + A+)V4, a, = 2(A+-1). (50) 

The value = 3 corresponds to the formation of a vacuum point at the trailing edge of the 
DSW [13] so that the density at the dark soliton minimum is Ug — = 0. For A'^ > 3 the 
vacuum point occurs inside the DSW at some x = Xy, where x~ < x^ < x~^ — see details in 



We define the relative intensity (hereafter - simply intensity) / of a DSW as the density 
ratio across it: 

(51) 

where ni and n2 are the values of density upstream and downstream the DSW respectively. 
This definition can be related to the one accepted in classical gas dynamics, where the 
relative pressure excess across the SW is often used as a measure of the SW strength. One 
should, however, stress that the notion of the DSW intensity for the NLS flows retains its 
original meaning only for DSWs not containing vacuum points. The modification of the flow 
resulting from the vacuum point appearance will be discussed below in Section 5.3.3. 

For the incident DSW (i.e. before the interaction) we obviously have rii = Us and ^2 = 1, 
i.e. its intensity is 

/«^^. (52, 

Now we turn to the left-propagating rarefaction wave, which is asymptotically described by 
the centred at x = / similarity solution of the classical limit equations (|5]), ([7]) (see, e.g., pT| . 
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(53) 



—— = V4\^,l) = , x,<x<xt; (54) 



Here the boundaries X2 are specified by the formulae: 



.Li 3-/i ~|~ X , , 

X2 = I — t , X2=l-\ 1 . (55) 



It is instructive to note that, since the modulation system fl20|) in the harmonic limit is 
consistent with the shallow- water equations ([5]), ([6]) — see fl27j) . the RW solution fl53l) . fl5^ 
is also a solution of full modulation system ( l20l) . namely 

A3 = A4 = A+; A2 = A+ = 1, Ai = A_(x,t). (56) 

This identification of the RW solution of the dispersionless limit equations as a particular 
solution of the full modulation system will be used in Section 5.2. The schematic behaviour of 
the Riemann invariants during the first stage of evolution, before the DSW-RW interaction, 
is shown in Fig. 7. 

One can readily see that dx^/dt > dx2/dt (this also follows from the characteristic 
velocity ordering described in Section 4.1) so the DSW will start overtaking the RW at some 
moment t = to when the leading edge of the DSW will catch up the trailing edge of the RW 
at Xo = xfito) = X2 (to)- Using fHOl) and (155!) we obtain 

- _ 2{A+f - 1 , 



5.2 Interaction, to < t < t* 

At t = to the DSW enters the RW region so that at t > to a nonlinear interaction zone 
confined to the interval forms (see Fig. 8) and evolves until some moment t* when 

the DSW completely overtakes the RW so that x^(t*) = x];{t*). At t > t* the DSW and 
RW fully separate, each acquiring a new set of parameters \j compared to their initial 
characterization. One should stress that, for t > to the functions xf{t) and xf{t) are no 
longer described by the formulae (H9|) . (!55l) from the previous subsection. 

The corresponding interaction diagram in the (x, t)-plane is shown in Fig. 9 (left). One 
can see that the NLS DSW-RW interaction in the semiclassical limit is essentially described 
by the interaction of two rarefaction fans: one of the shallow-water equations and another 
one — of the Whitham equations. 

In the interaction region [xg , x^] one still has A2 = 1 and A4 = but the remaining 
two Riemann invariants (Ai and A3) now vary so the modulation solution is no longer self- 
similar and a more general, hodograph solution ([31]) is needed. This is found via additional 
transformation (!36l) reducing Tsarev's equations (135!) for Wi^^i^Xi, X3) = Wi^^i^Xi,!, X3, A~^) 
to the EPD equation (1381) . The general solution (139|) of the EPD equation is parametrised by 
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Figure 8: Schematic behaviour of the Riemann invariants during interaction of the DSW 
and RW, to < t < t^. 




Figure 9: DSW-RW interaction diagram. Left: physical, {xt) plane; Right: hodograph, 
(A1A3) plane. 



two arbitrary functions 0i,2(A) which should be found from appropriate boundary conditions. 
These conditions, in their turn, must follow from the continuity matching conditions for Ai 
and A3 at the unknown boundaries X2{t) and x^it). 

At the left boundary x = X2 (t) of the interaction zone (segment PQ in the interaction 
diagram in Fig. 9, left) we have 

Ai = -1, A2 = l, A3 = A^(x2(t),t), A4 = A+, (58) 

where Xl{x,t) = X^lx/t) is found from the similarity modulation solution (H7|) . At the right 
boundary x = xf{t) of the interaction zone (segment PR in Fig. 9 left) we have, similar to 
the second condition (12^ . 

X^ = XL{xt{t),t), A2 = l, X, = X, = A+, (59) 

and X^_{x,t) = A_(^) is found from the rarefaction wave solution (15^ . 

We now need to translate nonlinear free-boundary conditions fl58p and (I59p into the 
boundary conditions for the function g{Xi, A3) satisfying the EPD equation 

2(A3 - Xi)dl,g = d,g - O^g , 9, = d/dXj . (60) 
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This is done in two steps. First we derive the boundary conditions for the functions 
A3) and W^i^Xi, A3) defining the hodograph solution 



x~Vit = Wi, x-Vst = W3. (61) 

Using the boundary condition at x = (!59|) and expression (1271) for the characteristic 
velocity Vi in the degenerate case when A3 = A4, the first equation fl6T]) becomes 

0\ _j_ 1 

x-^^t = W^i(Ai,A+). (62) 

Since according to the matching condition f l59|) one has Xi = X- at x = xf , we get, by 
comparing ( l62l) with the rarefaction wave solution ( [53l) . that 

iyi(Ai,A+) = /. (63) 

Next we turn to the boundary condition fl58l) and deduce from the comparison of second 
equation fl6T|) with similarity solution that 

W3(-1,A3) = 0. (64) 

Thus, the unknown at the onset curvilinear interaction zone PQTR in the (x,t)-plane maps 
to the prescribed rectangle PQTR in the hodograph (A1A3) plane (Fig. 9, right) exactly as it 
happens in the problem of the interaction of two simple waves in classical gas dynamics (see 
e.g. |4^). We also note that, in contrast to the original free-boundary matching conditions 
( 158|) . ( l59i) for the Riemann invariants Xj{x,t), the boundary conditions for the functions 
W^i,3(Ai, A3) are linear (i.e. they do not depend on the particular solution). 

To deduce boundary conditions for the EPD equation (!60!) from conditions (l63l) . (16^ for 
the Tsarev equations (135|) we use the relations (136|) between W^i,3(Ai, A3) and 5'(Ai, A3). Then 
from (163|) we obtain a simple ODE 

which is readily integrated to give the boundary value of the function g{Xi, A3) at A3 = A~^: 

^(Ai, A+) = Cii:(Ai, 1, A+, A+) + l= J^' + I , (66) 

where Ci is an arbitrary constant. 
Next, from ([MD, dSSD we find 

so the solution is readily found as 

(7(-l,A3) = C2£(-l,l,A3,A+), (68) 
where C2 is another arbitrary constant. 
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Conditions (166|) and (l68|l represent the Goursat type characteristic boundary conditions 
for the EPD equation (1601) . Now, we have two arbitrary functions 0i,2(-^) (see general solution 
(159]) ) and two arbitrary constants Ci_2 at our disposal to satisfy boundary conditions 
and flM|) . We first observe that, according to Section 4.4., the function g{—l,X3) has the 
meaning of the modulation phase shift in the incident DSW. Since this DSW is described 
by a centred simple wave modulation solution, this phase shift must be equal to zero. Thus 
we set C2 = so that condition (l68ll assumes the form 

g{-l,X,)=0 (69) 

in accordance with the phase shift requirement (H6|) . 

Now, the easiest way to proceed is to put 4>2{X) = and ai = — 1 in (15^ so that the 
solution of the EPD equation (J60l) reduces to a single quadrature 



(t>i{X)dX 



-1 



v/(A3-A)(Ai-A)' 



(70) 



Now we need to find 0i(A) and Ci to satisfy two conditions ( l66ll and ([69 
Substitution of (ITOll into boundary condition (J66i) yields 



W-*^ - ^' +/. (71) 



^(^+-A)(Ai-A) v'TF^ 
which is an Abel integral equation for 0i(A) (see e.g. [2]). We recall that 

if / -pS=d^ = fix) , then 0(x) = - / ^^^e 

a a 

Thus, the solution to ( 17T|) is readily obtained in the form 



Now one can see that condition fl^ is satisfied by flTU]) . fl72]) only if </)i(— 1) = 0, which 
implies that Ci = —ly/ + 1 and so finally 

Ai 



gf(Ai,A3) = [ — ^ dX 

' ' v/(^+-A)(A3-A)(Ai-A) 



TT 



(73) 



+ :(n,(.,.)-K(.)) 



7rv/(A+-Ai)(A3 + l) 

where ni(s,z) is the complete elliptic integral of the third kind (see, e.g. [2]) and 

^ (A+-A3)(Ai + l) Ai + 1 

^ (^+-Ai)(A3 + l)' ' A+-Ai- 



(74) 
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Hence, the modulation solution describing the interaction of counter-propagating DSW and 
RW is given by the formulae 



X, = 1, A4 = A+, X - Vi,3(Ai, 1, A3, A+)t = (^1 - ^^1.3) ^?(Ai, A3) , (75) 

where g{\i,\3) is specified by ( 173|) . 

The interaction continues until the moment t* defined by the condition a;^(t*) = (t*) 
(the right edge of the RW coincides with the trailing edge of the DSW - the point T on the 
diagram in Fig. 9, left). It is clear from the Riemann invariant sketch in Fig. 8 that this 
will take place when one has A3 = 1 and Ai = A~ simultaneously (see also Fig. 9, right). 
Substituting A3 = 1 and Ai = A~ into hodograph solution (1751) we find after some algebra 
that 

r= ^^"^"■) where . = i^llllf^ . (76) 

7i{l-A-)VA+-A-' 2{A+-A-) ^ ^ 

The corresponding coordinate x* = a;^(t*) = x~[{t*) is given by 

x*=(l + ^^—^] t* + P(l) , (77) 



2 

where -P(l) = g{A'', 1) (to be discussed below in Section 5.4). 
5.3 After interaction, t > t* 

At t = t* the DSW exits the RW region and the two waves separate. 




Figure 10: Schematic behaviour of the Riemann invariants after the interaction of the DSW 
and RW, t > U. 



5.3.1 Refracted DSW 

The modulation solution describing the DSW after the separation is given by three constant 
invariants (see Fig. 10) 

\i = A-, As = 1 , A4 = , (78) 
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while for the remaining one, A3, we have from (175|) a simple- wave modulation solution (cf. 

dSD) 



x = V3{A-,l,X3,A+)t + P{X3) 

(79) 



(A+-A3)(A3-1) 
(A3-l)-(A3-A-)/i(m) 



i(l + A- + A+ + A3) + ' ) t + ^(^3 



where 

(1-A-)M+-A3) , , 



and the function P(^) is found as 



7rv/(A+-A-)(e + l) V ' ' ' (e-A-)[(e-l)-(A+-l)Ml/)] 

(81) 

where 

A- + 1 A- + 1 A+-^ _ (l-A-)(A+-0 . . 

^ e + 1 ' ^ (A+-i)(e-A-)' ^ ^ 

Expressions (15^ are obtained from formulae (^^, where one sets Ai = A~, A3 = ^, and 
the modulus m in flHTl) is specified by flHOl) where A3 is replaced by ^. Thus, as a result of 
the interaction, the DSW is no longer described by the similarity modulation solution in 
the form of an expanding centred fan but rather becomes a general simple wave solution 
of the modulation system corresponding to the following initial-value problem for the NLS 
equation Q: 

A_(x,0) = A-, X+{x,0) = p-\x) , (83) 

P^^{x) being inverse of the function x = P(A+). The function -P(A+) in fl79p represents 
the DSW de-centring distribution acquired as a result of the interaction with the RW. It is 
directly related to the modulation phase shift distribution 6o{x,t) via fH5l) . fH6|) . It is not 
difficult to verify that P[^) = for = — 1. This is exactly what one should expect since 
when A" = — 1, there is no RW is generated and, therefore, there is no DSW refraction. 

The boundaries x^ and x^ of the refracted DSW are found by setting in fl79|l A3 = 1 (i.e. 
171 = 1) and A3 = A'^ (i.e. m = 0) respectively 

The background density and the amplitude of the trailing dark soliton in the refracted DSW 
are (cf. ([50])) 
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The intensity Ir of the refracted DSW is determined from (15T]) where we set rii = rig 
and n2 = |(1 — A~y (the latter is defined by the initial conditions (ITT]) ). Thus 



^r-[-Y^) ■ (86) 

5.3.2 Refracted RW 

The solution for the refracted RW is found from the hodograph modulation solution ( !75|) 
by setting in it A4 = A+ = A3 = A2 = 1, Ai = A_ (see ( l56ll ) and using that 

^i(Ai, A3, A3, A4) = V-{Xi, A4) (see EH]). As a result we get 

Q\ _|_ A + 

A+ = , x = l^_(A_, A+)t + ^(A^) = 1 + ^(A,) , (87) 

where the function G(^) has the form 

1V2 



(88) 



[(A+ + l)(ni(ri, r) - K(r)) + 2E(r 



where 



(A+-l)(f + l) f+1 

Similar to the refracted DSW, the refracted RW is no longer described by a centred fan 
solution but rather by a general simple- wave solution of the shallow- water system (I5]), (I6]) 
with the 'effective' initial conditions A+ = 1 and A_(x, 0) given by the function inverse to 
the refraction shift function G(A_). 

The boundaries of the refracted RW are given by the expressions 

X- = ^-^t + G{-1) , xt = - J t + G{A-). (90) 



5.3.3 Vacuum points 

As already was mentioned, an important property of the DSWs in the defocusing NLS flows 
is the possibility of the vacuum point (s) occurrence in the solutions for the problems not 
containing vacuum states in the initial data [H] . This effect has no analogue in both viscous 
SW dynamics and in the DSW dynamics in media with negative dispersion supporting 
bright solitons. Across the vacuum point, the flow speed changes its sign, which implies the 
generation of a counterfiow. As a result, the DSW with a vacuum point inside it, unlike 
a regular DSW, no longer represents a single oscillatory wave of compression: the vacuum 
point separates the compression part propagating to the right and the oscillatory rarefaction 
wave propagating to the left [H] . The DSW counterfiow due to the vacuum point occurrence 
has been recently observed in the experiments on nonlinear plane wave tunneling through 
a broad penetrable repulsive potential barrier (refractive index defect) in photorefractive 
crystals |56] . 
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If we fix the state ni = 1, Mi = in front of the DSW (as we do for the incident wave), 
then, by increasing the density jump n2 across the DSW we will be able to increase the 
DSW relative intensity only up to the value / = 4 at which the vacuum point occurs at 
the DSW trailing edge [23]. If n2 increases further, beyond the vacuum point threshold, 
the relative intensity of the compression part of the DSW decreases and, asymptotically as 
n^jnx —7- oo, vanishes so that the DSW completely transforms into the classical (smooth) 
left-propagating rarefaction wave [H]. This limit can alternatively be achieved by keeping 
the upstream state n2 fixed and letting n\ — )• 0: then we arrive at the well-known solution 
of the classical shallow- water dam-break problem (see e.g. [58]). 

Setting = — 1 we recover the already mentioned criterion > 3 for the vacuum point 
occurrence in the incident DSW. If which by (1851) . yields the relation A'^ = 2 — A , 

then the condition for the vacuum point appearance in the refracted DSW assumes the form 

A+ >2-A- . (91) 

The regions of the A~ , A'^ plane corresponding to different (with respect to the vacuum point 
appearance) flow configurations arising in the initial- value problem ([3]), ( JTOl) are presented in 
a diagram shown in Fig. 11. A particular flow evolution corresponding to Region II is shown 
in Fig. 12. We stress that, although the vacuum point appearance modifies the oscillatory 
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Figure 11: Regions in the plane of initial parameters {A~,A'^) — the classification with 
respect to the vacuum point occurrence. (I): No vacuum points; (II): No vacuum points 
in the incident DSW, a vacuum point in the refracted DSW; (III): Vacuum points in both 
incident and refracted DSWs. 

DSW profile (the lower DSW density envelope becomes nonmonotonous and the velocity 
profile acquires a singularity at the vacuum point — see [H], [31]), all the dependencies 
of the DSW edge speeds, density jumps and trailing soliton amplitudes on the initial data 
A'^,A~ remain unchanged. 

5.4 Key parameters of DSW refraction 

It is convenient to characterise the DSW refraction by three key parameters: the ampli- 
fication coefficient u which culd be defined as the ratio of the relative intensities f lSTj) of 
the refracted and the incident DSWs, the acceleration coefficient a which we define as the 
difference between the values of the DSW trailing dark soliton speeds s~ after and before 
the interaction, and the refraction shift d which is naturally defined as the phase shift of the 
DSW trailing soliton due to the DSW interaction with the RW (see Figs. 4,9). 
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Figure 12: Evolution of the profile ([TO]) with = 0, = 2.2, / = 50 (Region II in Fig. 11) 
leading to the occurrence of a vacuum point in the refracted DSW. 



For the first two parameters we readily have: 

2{A+ - A- 



'l-A-){l + A+) 



(92) 



see dHlD, ([52]), and 



a 



Sr — Sq 



dt 



dx^ 



t>t* 



dt 



l + A- 



> 0. 



(93) 



t<to 



- see iH), (|49]). 

Here the subscripts '0' and 'r' refer to the incident and refracted waves respectively. Note 
that the determination of v and a actually does not require knowledge of the full solution: 
both quantities are determined by the the transfer of the Riemann invariants through the 
DSW region. Interestingly, the acceleration coefficient a does not depend on the DSW 
strength before the interaction (~ A^) and is completely determined by the initial jump A~ 
of the Riemann invariant A_ across the RW. It also follows from ( l93l) that, since A~ > — 1, 
one has a > i.e. the DSW is always accelerated as a result of the head-on collision with 
the RW (indeed, a > implies acceleration of the trailing edge of the DSW and, therefore, 
acceleration of the DSW as a whole). The SW acceleration in the head-on collision with RW 
is also always the case in classical gas dynamics (see, e.g. |11]) as the SW meets the gas of 
decreasing density. 

Unlike the acceleration coefficient a, the amplification coefficient u can have both signs 
depending on the specific values of A'^ and A~ chosen, the boundary between the regions 
of the DSW (relative) strengthening and attenuation being given by equation A'^ = (1 — 
A^)/{1 + A^). We also note that, while the amplification coefficient u is formally defined 
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for the full range of values of A'^ and A , its original significance is retained only for the 
DSWs not containing vacuum points (see the discussion in the previous Section). 
The function (see (EI]), dHl) 

d{A\A-) = P(l) = ^-^^±^{U,{p,z*)-K{z*)) , (94) 

^{A+ - A-) 

where (see (l82l) ) 

A- + 1 A+ -1 A- + 1 

describes the refraction shift (see Fig. 4) of the trailing dark soliton in the DSW as a function 
of the initial parameters A~ . As a matter of fact, the determination of the refraction 
phase shift does require the knowledge of the full modulation solution in the interaction 
region. One can observe by comparing (IMl) with solution g{\i,\^) (1751) . ([71]) of the EPD 
equation for the DSW-RW interaction region, that 

d{A+,A-)=g{A-A), (96) 



which corresponds to the value of g at the moment t = t* (see f[76]) ). when the DSW exits the 
interaction region — see Figs. 8, 9. This is, of course, expected from the general modulation 
phase shift consideration described in Section 4.4. 

The dependencies of the refraction phase shift d on A" and A'^ given by ([94]), along 
with direct numerical simulations data for the refraction shift, are presented in Fig. 13. 
One can see that the dependence of the refraction shift on the density jump across the 
RW (roughly proportional to the value of A~) is much stronger than on the incident DSW 
strength (proportional to A'^). Along with the curves for the exact analytical solution for 
d, we also present the values of d extracted from the numerical simulations. One can see 
that, despite the fact that the accuracy of the modulation solution in the interaction region 
is not expected to be very high for the moderate spacing I between the initial jumps for the 
dispersionless Riemann invariants X±, the agreement between the asymptotic solution and 
direct numerics is quite good. The plots and comparisons with numerics for a and u will be 
presented in the next section as particular cases in the study of the DSW-RW interaction in 
the framework of a generalised, non-integrable version of the NLS equation. 




Figure 13: Typical behaviour of the DSW refraction phase shift d. Left: dependence d{A^) 
for fixed A'^ = 1.5, right: dependence d{A'^) for fixed A~ = 0. Solid line: formula ( |94]) . 
circles: direct numerical simulations data. 



One can trace certain analogy between the considered DSW-RW interaction and the 
two-soliton collisions in integrable systems: both interactions are elastic in the sense that 
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they both can be interpreted in terms of the "exchange" of spectral parameters by the 
interacting waves so that the global spectrum in the associated linear scattering problem 
remains unchanged. In the DSW-RW interaction the role of isospectrality is played by the 
transfer of the constant values of appropriate Riemann invariants of the modulation system 
through the varying DSW and RW regions so that one can predict the jumps of density 
and velocity across the refracted DSW and RW without constructing the full modulation 
solution. At the same time, the DSW and RW do not simply pass through each other 
and "exchange" the constant Riemann invariants: there are additional phase shifts for both 
interacting waves, similar to the classical soliton phase-shifts. The determination of these 
phase shifts requires knowledge of the full modulation solution. 



6 Refraction of dispersive shock waves in optical media 
with saturable nonlinearity 

6.1 Formulation of the problem 

We now consider the NLS equation with saturable nonlinearity (hereafter called the sNLS 
equation) 

where 7 > is the saturation parameter. This equation describes, in a certain approximation, 
the one-dimensional propagation of a plane stationary light beam through a photo-refractive 
crystal (see e.g. [21], [8]). One should note that in the nonlinear optics context the role 
of the time variable t is played by the spatial coordinate z along the beam propagation 
direction while x is the transversal coordinate. If the saturation effect is negligibly small 
(7|'?/'P ^ 1), then the sNLS equation ( 1971) reduces to the cubic NLS equation ([T]). The 
Madelung transformation ([21) with e = 1 maps equation ( 1971) to the dispersive hydrodynamics 
system (cf. ([3l)), 

nt + {nu)^ = 0, 
f n \ f nl n,A ^ (98) 
Vl + 7^A V8^' 

Here n has the meaning of the light beam intensity and u is the local value of the wave vector 
component transversal to the beam propagation direction. A detailed study of the periodic 
solutions to ( l98l) can be found in [13|. In particular, the linear dispersion relation for the 
waves of infinitesimally small amplitude propagating against the constant background flow 
with u = Uq, n = riQ has the form 



u = uoino,Uo,k) = kuo ± k^j f^^^^^y + ^ , (99) 

where u is the wave frequency and k is the wavenumber. 

In the dispersionless limit, system ( 19 8 p can be cast in the diagonal form with the Rie- 
mann invariants \± and characteristic velocities V± expressed in terms of the hydrodynamic 



28 



variables n and u as 



u 1 y/n 

X± = — ± —— arctan ^/Jn, V± = u± . (100) 

2 ^ ^ l + jn 



When 7 — )■ expressions (llOOp go over to the shallow- water relationships ([6]), ([7]) (note the 
different normalization for the dispersionless Riemann invariants compared to that used in 

Similar to (fTOj) . we specify the initial conditions for (!98|) in terms of two steps for the 
hydrodynamic Riemann invariants A-t 

f for X < 0, , / \ f — 7= arctan ^7 for x < I, 

= 1 arctan 77 for x > 0; ^-(^'0) = [ for x > /, 

(101) 

where > ^ arctan ^7, and — arctan ^7 < A < ^ arctan ^7. The special values of 
A+ for a; > and A_ for x < I are chosen such that initially the DSW and RW will propagate 
into an undisturbed "gas" (indeed, one can readily see that n = 1, n = in the middle 
region < a; < / (cf. ffTTj) ). 

Our numerical simulations of the evolution fl98|) . f llOip showed that, for a broad range 
of initial data parameters A^, the qualitative DSW refraction scenario is the same as in the 
cubic NLS case studied in previous sections. The quantitative characteristics of the DSW 
refraction, however, now depend not only on the initial conditions but also on the saturation 
parameter 7 entering the sNLS equation. This dependence was shown in to be quite 
strong for isolated photorefractive DSWs. We mention that knowledge of the effects of the 
photorefractive saturation on the parameters of a DSW is especially important in the context 
of an all-optical modelling of BEG dynamics (see [55|). Thus the DSW-RW interaction for 
the sNLS equation deserves a special study. 

Since the sNLS equation fl98l) is not integrable by the 1ST, the Riemann invariants are 
not available for the associated Whitham system and the modulation solution cannot be 
constructed by the methods used in previous Sections. An analytic description of the DSW 
refraction requires now a different technique. We shall take advantage of the theory of 
DSWs in photorefractive media developed in [I3] and based on the 'dispersive shock fitting' 
method introduced in [12] . As already was mentioned, our specific interest here is to quantify 
the effect of the nonlinear saturation on the DSW refraction, and, in particular, on the 
parameters a and u introduced above in the cubic NLS context (see (p3!) . (p2|) ). 



6.2 DSW transition relations 

The key ingredients of the dispersive shock fitting method of [12] in application to the sNLS 
equation (|98ll can be formulated as follows (see [13] for the details pertinent to the present 
study). Let the right-propagating DSW be confined to a finite region of space x~ < x < x~^ 
and connect two constant hydrodynamic states {ni,ui) ai x < x~ and (^2,^2) at a; > x~^] 
Hi > n2. Such a DSW is called a simple DSW. At the trailing edge x~ the simple DSW 
assumes the form of a dark soliton moving with constant velocity s~ and at the leading 
edge x~^ it degenerates into a vanishing amplitude linear wavepacket moving with constant 
group velocity s~^, > . The lines x^ = s^t represent free boundaries where the 
continuous matching of the mean flow (n, u) in the DSW region with the external constant 
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states (nijMi) and {n2,U2) occurs (in some cases it is more advantageous to formulate the 
matching conditions in terms of the mean density n and the mean momentum rm — see e.g. 



The simple DSW transition between the hydrodynamic states {ni,ui) and (^2,^2) is 
described by the following relationships: 



The value of the Riemann invariant A_ is conserved across the DSW , 



A. 



lx=x+ ' 



(102) 



I.e. 



(103) 



^-^ arctan Vl^, = ^ - — arctan ^ A° . 

The DSW edge speeds are defined by the kinematic conditions (cf. conditions ((30 
the cubic NLS case) 



on 

'dk 



n=n2, k=k'^ 



(104) 



n=ni, K=K 



The quantities k'^ (the leading edge wavenumber) and n~ (the trailing edge "soliton 
wavenumber" - the trailing soliton inverse half-width) in f ll04p represent the boundary 
values, k'^ = k{n2) and k~ = K{ni), of two functions k{n) and K,{n) satisfying the 
following ordinary differential equations (ODEs): 



dk 



dVt/dn 



dn v+{n) — dVt/dk'' 



dn 
dn 



dVt/dn 



v+{n) — dVt/dn 



fc(ni) =0; 



K{n2) = . 



Here 



v+{n) = V+{n,u{n)) = u{n) + 



n 



1 + 7n 



VL{n, k) = Uo{k, u{n), n) = k 
and 



uin) + 



n 



(1 + 7n)' 



+ 



A;2 



1 



uin) = 2 I A_ H arctan ^/yn . 

77 



(105) 
(106) 

(107) 

-iVL{n, IK,)] 
(108) 

(109) 



"Entropy" inequalities must hold ensuring that the hydrodynamic characteristics trans- 
fer data into the DSW region: 



Vl <s- < Vl, Vl < s- 



> s . 



;iio) 



Here V^. = V±{ni,ui), = V+{n2,U2) - see fllOOp for the definitions of V±{n,u). We 
note that inequalities f lllOp represent the dispersive-hydrodynamic analogs of classical 
Lax's entropy conditions [13]. 
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Relationships fll02p - flllOp enable one to 'fit' the DSW into the solution of the dispersionless 
limit equations without the knowledge of the detailed solution of the full dispersive system 
within the DSW region (much as in classical gas dynamics SW is fitted into the solution of 
the inviscid equations by means of the Rankine-Hugoniot conditions subject to Lax's entropy 
condition). 

Using the speed-amplitude relationship for the photorefractive dark solitons obtained in 
[T3] one can find the amplitude of the DSW trailing soliton. Setting the value s~ f ll04p 
of the DSW trailing edge for the soliton velocity c in formula (39) of [13] we obtain 



2{ni - as 



In 



1 + 7^1 



1 



,7«s 



l + l{ni- as) 1 + 7^1 



111) 



(note: Ui{ni) is given by the simple DSW transition condition (11030 ). 



6.3 DSW refraction 

Our concern in this section will be with the calculations of two DSW refraction parameters: 
the DSW amplification and acceleration coefficients, defined earlier in (p3|) and (p2|) as 

z/ = /.r//o and a = — Sq (112) 

respectively. We note that analytical determination of the refraction phase shift d is, un- 
fortunately, not feasible now as it requires knowledge of the full modulation solution, which 
is not available for the sNLS equation due to its nonintegrability so we shall present only 
numerical results for d. 



6.3.1 Before interaction, t < to 

The previous analysis of [13] suggests that the decay of two spaced initial discontinuities 
fllOip for the hydrodynamic Riemann invariants X± would result, similar to the cubic NLS 
case, in a combination of a right-propagating simple DSW centred at x = and a left- 
propagating simple RW centred at x = /. Indeed, the simple DSW transition condition 
f ll03p is satisfied by the initial step at x = 0, which implies a single DSW resolution of this 
step (provided the "entropy conditions" f lllOp are satisfied - see \13\ for the justification); 
similarly, the jump at x = / with constant Riemann invariant A+ across it asymptotically 
produces a single left-propagating RW (see Fig. 14). Indeed, our numerical simulations of 
the sNLS equation ( p8|) for a range of the saturation parameter 7 values confirm this scenario 
producing the plots qualitatively equivalent to that presented in Fig. 3. 

Now, following [L3l, we derive the key parameters of the simple photorefractive DSW in 
the form convenient for the further application to the refraction problem. 

To take advantage of formulae fll04p - fllOOp for the speeds of the DSW edges we first 
need to find the constant states {ni,ui) at x < x~ and (^2,^2) at x > x"*" defining the 
hydrodynamic jumps across the DSW. These are readily found from the the initial conditions 
fllOip and the relationship fll03p for the transfer of the Riemann invariant A_ across the 
simple DSW. According to the initial conditions fllOip the simple DSW must connect two 
hydrodynamic states with the same A_ = — j= arctan ^/7 while A+ = A'^ for x < x~ and 
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Figure 14: Distribution of the classical (dispersionless limit) Riemann invariants before the 
DSW-RW interaction 



A+ = ■^arctany7 for x > x+ (see Fig. 14). Then, using fllOSp and expressions fllOOp 
relating the Riemann invariants and the hydrodynamic variables n, u we find 

1 n 2 M+77 + arctany7\ +1 
no = 1, Mo = U , Til = — tan — , Mi = A arctan a/7 . 

(113) 

Thus, the /q of the incident DSW defined by f lSTj) is simply 

r 1, 2 A\^7 + arctany7\ 

Jo = - tan I ^- ^ ^ 1 . (114) 

Next, from (11031) we have A° = arctan ^/7, which by (I109p yields ^(r;.) = -^(arctan ^Jyn— 
arctan ^7) and so completely defines, via ffTOTi) . ffTOSll . ODEs ffT05l) . ffTnHll . 

As was shown in [13], it is convenient to introduce a new variable a instead of k using 
the substitution 

/ ^2(1 + 7n)2 ^ ^ 



which reduces ODE (I106p to the form 

da _ (1 + 5) [1 + 37n + 25(1 - 7fi)] 
(in ~ 2fi(l + 7ra)(l + 25) 



5(1) = 1. (116) 



The form flllGp has an advantage of being a separable ODE when 7 = 0, which makes it 
especially useful for the asymptotic analysis for small 7 . Once the function 5(n) is found, 
the velocity of the trailing soliton is determined by Eqs. f llU4p . f llOSp as 



^0 



(arctan ^7777 — arctan -v/7) H a{n\) , (117) 

1 + 7^1 

where n\ is given by Eq. (11130 . 

The amplitude of the trailing soliton is given by speed-amplitude relationship fillip . 
Using (imp . (I117P and the relationship u\ = (arctan ^7ni — arctan y^) following from 



103p one can derive the condition of the vacuum point occurrence at the DSW trailing edge 
see |l3j): 

5(ni) = (118) 
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Condition (IllSp yields, for a given value of the saturation parameter 7, the value of the 
initial density jump rii (and, therefore, of the parameter — see flll3p ) corresponding to 
the vacuum point appearance at the DSW trailing edge. Say, for 7 = 0.2 this value of is 
about 2.18 (cf. the critical value = 3 for 7 = 0) 

In conclusion of this Section we present an asymptotic expansion of Sq for small 7. First, 
to leading order we have from (11161) a separable ODE 

7 = 0: =-_^, a(l) = l, (119) 

an in 



which is readily integrated to give 



2 

a{n) = ^~l = 5o(n). (120) 
\/n 



We now introduce 

5 = 5o + 5i. (121) 
Substituting (11211) into (I116P and assuming 5i ~ 7 for 7 ^ 1 we obtain to first order 

dai 5i 4-3v^27 

-7— = -77^ + ^ ai(l) = 0, (122) 

dn 2n 4 - ^Jn \/n 

Eq. f ll22p is readily integrated to give 

27 



'n 



|^3(n-l) + 16(v^-l) + 64 In ^ ^ ^ 



(123) 



Now, substituting (I12ip . (I123P into (I117P and using expansion of n\ (11130 for small 7 we 
obtain to first order 

So = ^^-^ + ^ (1^ [(^^)' + ^^i^^f + - 245] + 128 In + 0(7') • (124) 

As one can see, expression (I124p agrees to leading order with the cubic NLS result ( 149|) 
for the trailing edge speed. We also notice that our perturbation approach formally breaks 
down for > 7 because of the logarithmic divergence in Eq. ( I124p as f 7 (we note that 
such values of A"^ correspond to very large density jumps {ni/n2 > 10) across the DSW — 
see p]). 

Formulae (11131) . (IllTp defining all the key parameters of the simple photorefractive DSW 
have been compared in [13] with direct numerical simulations data for a wide range of values 
of the saturation parameter 7 and a very good agreement was found. 



6.3.2 After interaction, t > t* 

Relations ( ]102p - ( lllOp describe a simple DSW transition between two constant states so 
they are not applicable to the varying transition in DSW-RW interaction region. However, 
one should still be able to use these relations for the determination of the key parameters 
of the refracted DSW when the interaction is over, provided no new hydrodynamic waves 
(DSWs or RWs) are generated and the output pattern consists only of a pair of the refracted 
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DSW and RW separated by a constant flow. If this is the case, one can say that the DSW- 
RW interaction is hydro dynamically (or semi-classically) "clean" (elastic). We stress that 
some zero-mean radiation due to non-integrability of the sNLS equation may be present but 
the latter does not affect the hydrodynamic transition conditions across the refracted DSW 
and RW. 

The notion of a hydrodynamically clean DSW-RW interaction can be elucidated by re- 
visiting the defocusing cubic NLS equation case considered in the previous sections. The 
"elasticity" of the DSW refraction for this case can be deduced from the following proper- 
ties of the initial-value problem for the defocusing NLS equation with a piecewise-constant 
initial datum ( ITTj) : (i) the asymptotic (e ^ 1) solution at any moment can only contain 
genus-zero (RW) or genus-one (single-phase DSW) regions, so in the semi-classical limit it 
can be globally described by the single-phase averaged NLS-Whitham equations - see [6]; 
(ii) the defocusing NLS-Whitham system is hyperbolic; (iii) the unique combination of the 
output (refracted) waves is determined by the transfer of the appropriate dispersionless limit 
Riemann invariants across the genus-one (DSW) and genus-zero (RW) regions. As a result, 
the "clean" DSW-RW interaction diagram on the x-t plane has the form shown in Fig. 9 
(left). 

We note that rigorous justifications of properties (i), (ii) (see 0) are based on the pres- 
ence of the integrable structure whereas property (iii) can be deduced using the classical 
method of characteristics for hyperbolic hydrodynamic type systems and does not require 
the availability of the Riemann invariants for the Whitham system and, hence, does not rely 
on integrability of the original equation (see |12j). 

Since the sNLS equation is not integrable we can only assume properties (i) and (ii) and 
then apply the transition relation fll03p (property (iii)) to the refracted DSW to determine 
the values n = rii and u = ui in the 'plateau' region between the refracted DSW and RW. 
Then our assumptions (i) and (ii) in the context of the sNLS equation could be indirectly 
justified a-posteriori by the comparison of the analytically obtained rii and ui with the 
corresponding density and velocity in the direct numerical solution of the IVP flUS]) . fllOll) . 

Since the refracted DSW propagates to the right, into the region with A_ = A~ (see the 
initial conditions (110 ip at x — +oo) one must have, by (11020 . the same A_ = A~ across it, in 
the constant 'plateau' region. Next, the refracted RW propagates to the left, into the region 
with A+ = A'^ (again, see initial conditions (110 ip at x — )■ —oo) and, therefore, A+ = 
everywhere through this wave and in the 'plateau' region. From the initial condition f ll03p . 
the value of A_ to the left of the RW is A_ = — ^arctan.^ and the value of A+ to the 

right of the DSW is A_|_ = ^ arctan y/^. Thus, we arrive at the Riemann invariant diagram 
schematically shown in Fig. 15 (cf. diagram in Fig. 10 for the cubic NLS case). 

Thus, using relationships ( llOOp between the Riemann invariants X± and the hydrody- 
namic variables n,u, one arrives at the set of equations determining the hydrodynamic 
states {ni,ui) and (^2,^2) at the trailing and leading DSW edges respectively: 

Ui 1 , Ui 1 U2 1 

\ arctan A/7r?,i = A : arctan \hn\ = arctan ^ryno = A : 

2 ^ 2^7 ^' 2^7 ^ ' 

U2 1 , 1 _ 

1 arctan xrrru = arctan . 

(125) 
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Figure 15: Distribution of the dispersionless limit Riemann invariants after the DSW-RW 
interaction. 



So 



ni = -tan"fy^ 1, Ui = + A' 

2/1 . ^ A- ^\ . 1 



77-2 = — tan - arctan , U2 = A -\ arctan ^/y . 



(126) 



To verify our key assumption about the "semi-classically clean" DSW-RW interaction 
in the sNLS equation case we have compared the values of the density and velocity in the 
region between the refracted DSW and RW obtained from direct numerical simulations of the 
sNLS equation with the predictions for rii and ui of formulae f ll26p based on this assumption. 
As one can see from Fig. 16 the the comparisons show an excellent agreement confirming 
our clean interaction hypothesis for a range of values of 7, A'^ and A~ . At the same, one 
can notice a small discrepancy visible at larger values of A'^ > 1.7) in the plots for 
z/(A+). This is connected with the occurrence of the vacuum point in the refracted DSW 
for sufficiently large density jumps across it. As was observed in [13], for large-amplitude 
photorefractive DSWs the Riemann invariant transition condition f ll03p is replaced by the 
classical Rankine-Hugoniot shock jump conditions so relation fll26p holds for large A'^ only 
approximately. 




Figure 16: Density rii in the constant flow region between the refracted DSW and RW. 
Left: ni{A~) for fixed A'^ = 1.5; Right: ni{A'^) for fixed A~ = 0. Solid lines: analytic 
(modulation theory) curves; dots: direct numerical simulations data. 



Now, we shall use general relationships fll04p — fillip to derive the trailing soliton pa- 
rameters in the refracted DSW. 
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Comparing (11030 and (I125P we find A° = A so expression (I109P for u{n) assumes tlie 
form 

M(ra) = 2(y4~H arctanA/71) . (127) 



Substituting (11270 into (11070 and (11080 and using tlie same cliange of variable (11150 in 
ODE ( I1O60 we arrive at tfie same ODE ( 11160 for tfie function a{n) but now with a general 
boundary condition 5(n2) = 1 since ^2 7^ 1 for the refracted wave (see ( 11260 ). As before, 
this condition follows from the boundary condition for k in ( I1O60 and the relationship (IllSp 
between a and k. The velocity of the trailing soliton in the refracted DSW is determined by 
Eqs. (fTOill . (fT08|) as 



2 A 



arctan -v/T^ 1 + 

^ y 1 + 7ni 



(128) 



where rii is now given by Eq. (I126p . Comparison for the dependence s~{A~^) for a fixed 
value of = —0.8 is presented in Fig. 18. One can see that the value of s~ quite strongly 
depends on the saturation parameter 7. Expanding s~ for small 7 we get (cf. ( 1124p ) 




Figure 17: The refracted DSW trailing edge speed s~ as a function of an input parameter 
A'^ for fixed A~ = —0.8. Solid lines: modulation solution ( I128P : dots: numerical simulations 
data. 
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(129) 

2 , " — 2 ■ Again, one can see that the leading order of expansion (I129p 
agrees with the cubic NLS result ( 18^ as expected. 

Given the value of s~, the trailing dark soliton amplitude in the refracted DSW is 
found from formula (11110 . Comparisons of the analytically found values of as for 7 = 0.2 
with direct sNLS numerical simulation data are presented in Fig. 17. and show excellent 
agreement. Also, the dashed lines show the dependencies as{A~) and as{A'^) for 7 = 0. 
As one can see, the nonlinearity saturation has strong effect on the refracted DSW soliton 
amphtude. The condition = rii defining the vacuum point occurrence at the trailing 
edge of the refracted DSW, leads to the same equation ( 11180 . which was obtained earlier for 
the incident DSW, with the only (essential) difference that rii is now given by ( 11260 . The 
vacuum point regions diagram for 7 = 0.2 is presented in Fig. 19. 
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Figure 18: Trailing soliton amplitude a^. Left: as{A~) for A'^ = 1.5; Right: as{A'^) for 
A~ = —0.4. Solid line: analytic curve for 7 = 0.2; Dots: direct numerical simulations data 
for 7 = 0.2. Dashed line: the curve for 7 = 0. 




Figure 19: Regions of the plane of initial parameters A'jA'^ for 7 = 0.2: (I) No vacuum 
points; (II) No vacuum points in the incident DSW, a vacuum point in the refracted DSW; 
(III) Vacuum points in both incident and refracted DSWs. 
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Comparison with the analogous diagram for the Kerr nonhnearity case 7 = (Fig. 11) 
shows that variations of the saturation parameter 7 have rather significant effect on the vac- 
uum point appearance. Our numerical simulations confirm this conclusion. As already was 
mentioned, in the developed modulation theory we assume a semiclassically "clean" DSW- 
RW interaction, which, strictly speaking, applies only to the region I in Fig. 19. However, 
our comparisons show that, if the initial parameter A'^ is not too large, the DSW fitting ap- 
proach [12] implementing the Riemann invariant transition condition f llOSp gives reasonably 
good quantitative predictions for the refracted DSW parameters for the regions II and III 
as well. 



6.4 DSW refraction parameters 

The DSW amplification coefficient is defined as u = Ir/Io-, where the incident DSW relative 
intensity Jq is given by (11141) . Using (I126p the relative intensity of the refracted DSW is 
readily found in terms of the input parameters and A~ as (see fl5T|) ) 

tan2(y7^il^) 

Ir = — = ^71 ^ Tz r • (130) 

n2 taxr (^ arctan ^7 — ^y/l) 

In Fig. 20 we present the dependencies z^(v4~) (for a fixed A^) and z^(v4+) (for a fixed 
A~). One can see that the amplification coefficient (unlike the individual parameters of the 
incident and refracted DSWs — see e.g. Fig. 16 above and Figs. 18, 19 below) shows a very 
weak dependence on the saturation parameter 7 for rather broad intervals of A^ and A~ 
so that one can safely use simple expression fl92|) obtained for 7 = 0. The direct numerical 
simulations fully confirm this conclusion (we do not present numerical points on Fig. 20 to 
avoid cluttering the plot). 




Figure 20: DSW amplification coefficient v. Left: iy{A~) at = 1.5, v4+ = 1.5. Right: 
z/(A+) at A- = 0; 



Now we look at the behaviour of the acceleration coefficient a = s~ — Sq , which is found 
analytically with the aid of formulae fll28p and (11170 . The dependence (7(7) for A^ = 1.2, 
A^ = —0.7 (Region I in Fig. 19) is shown in Fig. 21. One can see that, similar to the 
amplification coefficient u, the dependence of a on 7 and A^ (i.e. on the intensity of the 
incident DSW) is quite weak. Indeed, the relative change of a does not exceed 10% over the 
broad interval of 7 from to 0.5). Thus, at least in region I, one can safely assume the simple 
expression (p3|) cr = (1 + A^)/2 obtained for the cubic nonlinearity case. The comparisons 
with numerics presented in Fig. 22 confirm this observation. To analytically quantify the 
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Figure 21: Analytical curve for the DSW acceleration coefficient cr as a function of the 
saturation parameter 7 for A'^ = 1.2, = —0.7. 
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Figure 22: The DSW acceleration coefficient a as a function of input parameters A and 
A^. Dashed lines: analytic curves for 7 = 0; Cirles: numerical data for 7 = 0.2. Left: cr[A~) 
at fixed A+ = 1.2; Right: a{A+) for fixed A' = -0.4 



deviations of the quite complicated general "photorefractive" dependence a{A'^ , A~ ,'-/) from 
the simple dependence a = {1 + A~)/2 in the cubic nonlinearity case given by f l93|) . we derive 
an asymptotic expansion for a for the case when both interacting waves have small intensity. 
Introducing and £_ by 

1 1 

A = arctany7 + e_, A'^ = arctany7 + e+ (131) 

and assuming e_ <^ 1, e+ ^ 1 we obtain from ( 1124^ and ( I129P on retaining second order 
terms, 

a = - So = Y + £-7 + 0(e_7^; e^r, (132) 

One can see that expansion (11321) does not contain terms proportional to e+7, which implies 
that, for the interactions involving weak photorefractive DSW and RW, the acceleration a 
of the DSW up to second order does not depend on its initial intensity. 

Finally, in Fig. 23 we present numerical values for the DSW refraction shift d (see Fig. 4) 
taken for the particular value of 7 = 0.3. The numerics (circles) are put against the analjd;ical 
curves d{A~,A~^) defined by formula fl94|) for the cubic nonlinearity case, 7 = 0. One can 
see that, similar to other definitive DSW refraction parameters u and a, there is almost no 
dependence on A'^ and 7 at a fixed value of A~ (roughly, the RW intensity), however, the 
departure of the dependence d on A~ from the Kerr case 7 = becomes more pronounced 
with growth of A^ . 
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Figure 23: DSW refraction phase shift d. Left: dependence d on for fixed A'^ = 1.5; 
Right: dependence d on for fixed A~ = 0. Dashed hnes correspond to 7 = (analytical), 
circles — to 7 = 0.3 (numerical). 

7 Conclusions 

In this paper, we have considered a dispersive counterpart of the classical gas dynamics 
problem of the interaction of a shock wave with a counter-propagating simple rarefaction 
wave often referred to as the shock wave refraction problem. Apart from the obvious contrast 
between both local and global structures of viscous SWs and DSWs, there is a fundamental 
difference between the classical dissipative, and the present, dispersive conservative settings 
which makes possible full quantitative description of the DSW refraction. The salient feature 
of the viscous SW refraction is the generation of the varying entropy wave resulting in 
a complicated system of the Rankine-Hugoniot shock conditions resolvable in most cases 
only by numerical means. Contrastingly, in conservative dispersive hydrodynamics, the 
thermodynamic entropy does not change and the jumps of the hydrodynamic quantities 
across the DSW are completely determined by the transfer of the Riemann invariants of 
the appropriate modulation Whitham equations along the characteristics. Essentially, the 
DSW-RW interaction problem reduces, in the semi-classical limit, to the description of the 
interaction of two expansion fans: one of the shallow-water equations and another one - of 
the Whitham modulation equations. 

Our study was performed in the frameworks of the one-dimensional defocusing NLS equa- 
tions with cubic nonlinearity (Eq. [T]) and saturable nonlinearity (Eq. W7}i . To model a generic 
DSW-RW bidirectional interaction we have considered the initial-value problems for both 
NLS equations with the initial data given by appropriate piecewise-constant distributions for 
the density (the wavefunction squared modulus) and the velocity (the wavefunction phase 
gradient). To single out the "pure" DSW-RW interaction, we specified the initial data in the 
form of two steps for the "Eulerian" (dispersionless limit) Riemann invariants having jumps 
of different polarity shifted with respect to one another by a large distance / (see Fig 2a). 

For the integrable cubic nonlinearity case we have constructed exact modulation solutions, 
asymptotically (t ^ 1) describing all stages of the bidirectional DSW-RW interaction in 
terms of the evolution of the Riemann invariants of the NLS- Whitham system. This was done 
by mapping the original nonlinear Gurevich-Pitaevskii type matching modulation problem 
to the Goursat problem for the classical linear Euler-Poisson-Darboux equation f l60p . Along 
with the modulation solution describing slow variations of the amplitude, the wavelength, 
the mean etc. in the DSW, we have derived explicit compact expressions for the DSW- 
RW refraction phase shifts, having certain analogy with the classical soliton phase-shifts in 
two-soliton collisions. 
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For the NLS equation with saturable nonhnearity, which is a typical model for the de- 
scription of the light beam propagation through photorefractive optical materials, we have 
taken advantage of the DSW fitting method [12j applicable to non-integrable dispersive sys- 
tems. This method was applied recently in [13j to the description of the simple-wave optical 
photorefractive DSWs and in the present study we extended it to the DSW-RW interac- 
tion. Our consideration of "non-integrable" DSW refraction in the framework of the NLS 
equation with saturable nonlinearity ( 1971) is based on the assumption (confirmed by direct 
numerical simulations) that the head-on DSW-RW interaction is "semiclassically elastic", 
i.e. is not accompanied by the generation of new DSWs or/and RWs. The comparisons of 
the key parameters of the photorefractive DSW refraction: the amplification coefficient v 
and the acceleration coefficient o defined by formulae (11121 a) and (11121 b) respectively, with 
their Kerr (7 = 0) counterparts have revealed a rather weak dependence of these particular 
parameters on the saturation coefficient 7, which could prove useful for the experimental 
all-optical modelling of the BEG DSW refraction using photorefractive materials. 

The direct numerical simulations of the photorefractive DSW refraction confirm key 
predictions of our modulation analysis, which provides further striking evidence of the ro- 
bustness of the modulation theory in non-integrable dispersive wave problems, now in the 
more complicated setting involving DSW-RW interactions. 

We conclude with the remark that the approach used in this paper can also be applied to 
obtain asymptotic solution to the problem of the overto/^m^ DSW-RW interaction in the NLS 
flows. While this problem was studied in the KdV equation framework in [1] , we believe that 
it deserves special attention in the context of the defocusing NLS equation since, due to a 
different dispersion sign and the possibility of the vacuum point occurrence within the DSW 
one can expect a number of qualitative and quantitative differences compared to the KdV 
flows. Also, the developed theory can be readily extended to the problem of the generation 
of DSWs by the nonlinear dispersive interference of two simple rarefaction waves studied 
numerically in |33] . 
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